Comparison of phase-field models for surface diffusion 
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The description of surface-diffusion controlled dynamics via the phase-field method is less trivial 
than it appears at first sight. A seemingly straightforward approach from the literature is shown 
to fail to produce the correct asymptotics, albeit in a subtle manner. Two models are constructed 
that approximate known sharp-interface equations without adding undesired constraints. Linear 
stability of a planar interface is investigated for the resulting phase-field equations and shown to 
reduce to the desired limit. Finally, numerical simulations of the standard and a more sophisticated 
model from the literature as well as of our two new models are performed to assess the relative 
merits of each approach. The results suggest superior performance of the new models in at least 
some situations. 
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I. INTRODUCTION 

For a large class of pattern-forming systems, the essen- 
tial dynamics to be understood and described is that of 
an interface between two phases. Mathematically speak- 
ing, part of the problem to be solved consists in deter- 
mining the position of the interface as a function of time, 
i.e., is a free or moving-boundary problem. 

Phase-field models have been established as powerful 
tools for the numerical simulation of this kind of prob- 
lem. They avoid explicit front tracking and are versa- 
tile enough to deal with topological changes. Phase-field 
methods constitute a special case of level-set approaches 
[TJ |2]. They differ from the general case by having a 
level-set function that satisfies particular bulk equations 
of motion rendering unnecessary the computation of an 
extension of the interface velocity to the bulk. This way, 
they also avoid interface capturing at each time step, 
which is normally requisite in level-set methods. In com- 
parison with other level-set methods, a drawback of the 
phase-field approach is that it does not yield an exact rep- 
resentation of the interface continuum problem, reducing 
to its dynamics only asymptotically. However, quantita- 
tive numerical control of phase-field representation has 
made enormous progress in the last decade [HH], so this 
disadvantage is not crucial anymore. 

In a phase- field model, information on the interface po- 
sition is present implicitly, given either as a level set of 
a particular value of the phase field (in two-phase mod- 
els) or by equality of the phase-field values for different 
phases (in multi-phase models) , and can be recovered by 



computation of the appropriate level set at only those 
times when knowledge of the position is desired. 

A major field of application are solidification problems, 
where diffuse-interface models were developed early on 
El El 13 and have seen renewed interest ever since com- 
putational power increased enough to render their simu- 
lation feasible. The concept was extended to anisotropic 
interface properties [5], and first qualitative numerical 
calculations of dendritic growth El [TU] were followed by 
theoretical improvement of the asymptotics permitting 
quantitative simulations [3JH], at least for intermediate 
to large undercoolings. Non-dendritic growth morpholo- 
gies were also simulated, even in three dimensions 11J. 
Generalizations included the description of the coexis- 
tence of more than two phases [TJ] [T3] . 

Additional examples of successful application of the 
tool phase field include the modeling of step-flow growth 
|14[ IT5] and of the elastically induced morphological in- 
stability HSHnjUH], often labeled Grinfeld QJ] or Asaro- 
TiUer-Grinfeld (ATG) instability [20]. All of the exam- 
ples mentioned so far dealt with nonconservative interface 
dynamics, where a particle reservoir is provided by either 
the melt that is in contact with the solid or the adatom 
phase on a vicinal surface. 

Actually, regarding the ATG instability, which is an 
instability with respect to material transport driven by 
elastic energy, interest initially focused on transport by 
surface diffusion, which leads to conserved dynamics. 
This was the case in the first article by Asaro and 
Tiller [20] . but also in the first numerical simulations by 
sharp-interface continuum models (21 j . preceding com- 
putations of the instability under transport by melting- 
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crystallization [22 . 

The situation reversed when the phase-field method 
was employed for the first time to compute the ATG in- 
stability [16j [17] . Here, all the early works considered a 
nonconserved phase-field [THl HZ] (T5J [23J . Only recently 
has surface diffusion been considered in phase-field mod- 
els treating elastically stressed materials [2H US] ■ This 
difference in preferences when modeling either on the ba- 
sis of a sharp-interface model or using a phase field may 
be due to the fact that writing down a nonconservative 
and a conservative model is equally simple in the former 
case, whereas it is less straightforward to write down the 
conservative model within the phase-field approach than 
the nonconservative one. 

This is not to say that phase-field models with a con- 
servation law for the phase field have not been considered 
at all. Starting from a Cahn-Hilliard equation with a con- 
centration dependent mobility, Cahn et al. |26j obtained 
an interface equation with the normal velocity propor- 
tional to the Laplacian of the mean curvature, ft then 
appears as if all phase-field models with surface diffu- 
sion should be derivable on the basis of similar consid- 
erations. Indeed, comparable models have been applied 
in the simulation of clcctromigration and voiding in thin 
metal films [27] [25]. These two models are slightly dif- 
ferent, but the difference is not crucial and all previous 
models except the one given in 24] seem to suffer from 
the same problem, to be discussed in the following. 

As we shall see, it is quite easy to set up a conservative 
phase-field model. But it is more difficult to obtain the 
correct asymptotics describing surface diffusion as given 
by the desired sharp-interface limit. Past models such 
as the ones presented in [25] [26] [2JJ [28] while asymp- 
totically producing a set of equations containing the de- 
sired limit equations, include an additional restriction, 
i.e., they have one equation too many, a fact that seems 
to have been overlooked so far. In [24], this restriction is 
not present, but the authors consider their improvement 
only a stabilizing element, not changing the asymptotics, 
whereas what they have achieved in reality is superior 
asymptotic behavior. Because the flaw of the faulty mod- 
els is subtle, it is not a priori clear how adversely the 
undesired restriction will affect their behavior. There- 
fore, numerical simulations are necessary to assess their 
respective virtues and drawbacks. 

The purpose of the present article is to demonstrate the 
overlooked restriction, to explore an alternative approach 
to phase-field modeling of surface diffusion, to derive ad- 
ditional models not having the aforementioned flaw, and 
finally, to compare the different models numerically. 

To render things as simple as possible, we will restrict 
ourselves to two dimensions and give analytic derivations 
only for purely surface-diffusion-driven motion, i.e., the 
coupling to a destabilizing process such as clastic relax- 
ation or electromigration will not be considered in the 
asymptotic analysis. The fully three-dimensional model 
including elastic energy and thus describing the ATG 
instability has been given in |29j . an article with lower 



pedagogical ambitions than this one. In simulations, we 
will consider both surface diffusion and elastic degrees 
of freedom, i.e., the ATG instability, to be able to make 
comparisons for stable and unstable situations. 

The paper is organized as follows. In Sec. [TT] the 
sharp-interface model to be approximated by the phase- 
field equations is specified. The nonconservative case 
will be discussed for reference purposes. Section [f f f | then 
presents the standard approach that previously was sup- 
posed to reduce to the correct limit and pinpoints the 
oversight in existing asymptotic analyses. An alternative 
approach is presented in Sec. |1V| failing for complemen- 
tary reasons. By appropriate combination of the ideas 
from both approaches, two phase-field models will be 
given in Sec. \V\ producing the correct asymptotic behav- 
ior. An analytic linear stability analysis of these models 



shows, in Sec. VI that they correctly reproduce the spec- 
tra of the sharp-interface limit. In Sec. |Vff| comparisons 
of the different models will be performed via numerical 
simulation for a number of pertinent situations. Finally, 
some conclusions to be drawn from both analytic and 
numerical calculations will be discussed. 



II. SHARP-INTERFACE MODEL FOR MOTION 
INDUCED BY CURVATURE 

In the simplest case, where surface energy is the only 
relevant quantity determining the motion of an interface, 
the local chemical potential difference between the solid 
and the second phase [liquid, gas (vapour), or vacuum] 
at the interface may be written 



1 

8fi= —jk, 

Ps 



(1) 



where p s is the density of the solid phase, y the (isotropic) 
surface tension and k the curvature (in 3D, the mean 
curvature) . A positive curvature corresponds to a locally 
convex solid phase, a negative one to a locally concave 
solid. 

Once the chemical potential is given, the stability of 
the interface can be assessed. From ([lj, we infer immedi- 
ately that a planar interface is energetically stable. Any 
protrusion of the solid gives a convex bump and increases 
the energy of the solid over that of the second phase, lead- 
ing to a noncquilibrium situation favoring diminution of 
the solid phase. Any indentation of the solid produces a 
concave bump and decreases the energy of the solid be- 
low that of the second phase, leading to a noncquilibrium 
situation favoring growth of the solid phase. 

However, in order to determine the evolution of an un- 
stable state, some dynamical law governing its motion 
must be stated. If a particle reservoir is present and the 
interface is rough, it is natural to assume linear noncqui- 
librium kinetics. The driving force then is the chemical 
potential difference itself, and the normal velocity v n of 
the interface will be proportional to it: 



v„ = -k v 3fi , 



(2) 



3 



where k v is a mobility and the normal points from the 
solid into the second phase. 

On the other hand, for material transport by surface 
diffusion, the driving force is the gradient of the chemical 
potential along the surface, producing a surface current 
j cc -V s 6fi (V s is the surface gradient), which leads to a 
dynamical law of the form 



M s A.,<5/i = M s 



d 2 5n 
"ds 2 



(3) 



where A s is the Laplace-Beltrami operator on the surface, 
reducing to a double derivative with respect to the arc 
length for a one-dimensional interface, and M s a mobility 
coefficient (dimensionally different from the mobility k v ) , 
assumed constant here. 

A linear stability analysis of a planar interface is read- 
ily performed, writing 



ikx+iot 



(4) 



where C,q is the constant position of the unperturbed in- 
terface, x the cartesian coordinate parallel to it, and e 
a small parameter used to keep track of orders of the 
perturbation expansion. The form of the perturbation, 
containing a wave number k and a growth rate a>, is dic- 
tated by the fact that plane waves are eigenmodes of the 
differential operator d 2 appearing in the definition of the 
curvature ^ and that v„ = £/[l + (d x ^) 2 ] ] ^ 2 is propor- 
tional to a first-order time derivative. Using 



[1+(<U) 2 ] 3/2 ' 
one obtains the dispersion relations 

Ps 

for the nonconservative and 

a) = k = -Mk 



(5) 



(6) 



(7) 



for the conservative cases, respectively. To arrive at the 
last result, note that to linear order (in e) d 2 s is not dif- 
ferent from d 2 . 

For brevity, we have defined new kinetic coefficients K 
and M, which allows us to avoid carrying along the factor 
ylp s all the time. 

The two models for motion by curvature considered 
here are given by Eqs. ([!]) and ^ on the one hand and 
Eqs. ([TJ and ^ on the other, describing the noncon- 
servative and conservative cases, respectively. A phase- 
field model trying to represent these dynamics should 
converge to the appropriate set of these sharp-interface 
equations in the limit of asymptotically small interface 
width. More interesting and more complex physical prob- 
lems are obtained, when the normal velocity depends on 
additional ingredients beyond curvature-induced driving 
forces. This will then lead to a coupling of the phase field 
to other fields. We shall discuss such a case later on. 



III. SCALAR-MOBILITY PHASE-FIELD 
MODEL 

Before considering the structure of previous phase-field 
models attempting to capture surface diffusion dynamics, 
let us briefly recall the phase-field model for nonconserved 
order parameter (p. This can be written |18j 



dt 



(8) 



where f((f>) - (p 2 (l - <p) 2 is the usual double-well potential 
describing two-phase equilibrium. cf> varies between 0, 
corresponding to the nonsolid phase, and 1, correspond- 
ing to the solid; W measures the width of the transition 
region between the two phases, i.e., it may be interpreted 
as an interface thickness. A prime denotes a derivative 
with respect to the argument. 

The standard approach to a phase-field description of 
surface diffusion, as proposed in [25J ED EZ1 12H], is then 
to prepend the right hand side of Eq. (8j with a dif- 
ferential operator corresponding to the divergence of a 
gradient multiplied by a phase-field dependent mobility, 
i.e., Eq. (|8| becomes replaced with 



j = MV^<5/i(V 2 <^) 
#t(V 2 cf>, 4>) ee - W 2 V 2 (p + 2f(<p) 



(9) 



where M is a scalar function of either <p (25J HS1 HH] or 
WV<p [27 , chosen such that the mobility tends to zero far 
from the interface: M(<p, WV(f>) — > for <p — > and <p — > 1. 
5jx is a nondimensionalized chemical potential difference. 

We will refer to the model described by Eqs. ^ as the 
scalar-mobility model or briefly SM model in the follow- 
ing. 

At this point, a few remarks are in order. First, the 
field cp is the density of a conserved quantity by con- 
struction, since the right hand side of ([9]) is written as a 
divergence. This is true for any (nonsingular) form of the 
mobility function M. Second, 5jl becomes zero for <p — > 
and — » 1 , meaning that there is no diffusion in the bulk 
anyway. One might therefore wonder whether it is really 
necessary to choose a mobility that goes to zero in the 
bulk. The conservation law plus the absence of diffusion 
far from the interface should suffice to restrict transport 
to diffusion along the interface. In fact, we shall see that 
essentially the same asymptotic results are obtained no 
matter what the form of M, the only conditions to be im- 
posed being positivity (for almost all values of (p or V<p) 
and boundedness. It is just easier to derive them if it is 
in addition assumed that M vanishes in the bulk. On the 
other hand, it will turn out that if a restriction imposed 
by the asymptotics is removed (or not yet satisfied in the 
temporal evolution of the system), M has to decay suf- 
ficiently fast inside the bulk for the limit to make sense. 
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This may be relevant for the behavior of the model before 
it reaches its asymptotic state. 

Finally, the issue at present is not so much whether 
the dynamics is conservative but whether it does reduce 
to the sharp-interface model of Sec. [IT] in the limit of an 
asymptotically vanishing interface thickness. To investi- 
gate this, we have to explicitly carry out the asymptotic 
analysis. 



& a &p, where a, (3 e [r, s). The metric tensor reads 



{gap) - g - ( n n 



(1 +rK) z 



(12) 



its determinant is 



g = detg = (1 + yk) 2 , 



(13) 



A. Local coordinate system 

The basic idea of the analysis is to expand all dynam- 
ical quantitities in terms of the small parameter W de- 
scribing the interface thickness, to solve for the phase 
field and to use the solution to eliminate its explicit ap- 
pearance from the equations. To this end, the domain 
of definition of the field is divided into an outer region, 
where gradients of fields can be considered to be of order 
one and an inner region (close to the interface), where 
these gradients are of order l/W. The expansion in pow- 
ers of W is rather straightforward in the outer domain, 
Eq. @ can be taken as a starting point directly. As to the 
inner domain (and its matching with the outer region), it 
is useful to first transform to coordinates adapted to its 
geometry. Therefore, local coordinates r and s are intro- 
duced with r orthogonal to the interface (defined as the 
level set corresponding to <p(x,z,t) — 1/2) and s tangen- 
tial to it. r is the signed distance from the interface and 
will be rcscalcd by a stretching transformation r = Wp 
to make explicit the W dependence for the expansion, s 
is the arc length of the interface curve. Inner and outer 
solutions must satisfy certain matching conditions due 
to the requirement that they agree in the combined limit 
W — ♦ 0, p — > ±oo, r — * 0. These conditions are given in 
App.0 

To obtain a set of basis vectors for our local coordinate 
system, we first write 



r = R(s) + rn(s) , 



(10) 



where r is the position vector of a point near the in- 
terface, R the position of the interface itself, and n the 
normal vector on it (oriented the same way as in the 
sharp-interface model, i.e., pointing out of the solid). 

Given the coordinates, it is a trivial matter to write 
down a coordinate basis 



or 

dr <9R <9n „ 

fi * = IT = IT +r lT =a + «)t. 
as as os 



(11) 



which is orthogonal. (This is no longer automatically 
true in 3D |29j.) t = dR/ds is the unit tangent vector to 
the interface, and dn/ds = Kt is one of the Frenet formulas 
[30J , specialized to two dimensions. 



hence yfg - 1 + rn (using the locality of r to ascertain 
tk < 1 ) , and the contravariant components of the metric 
tensor are obtained as 



{g af) ) g 1 | (i + «r 2 )' 



(14) 



From now on, we use the Einstein summation con- 
vention for pairs of covariant and contravariant indices. 
The vectors of the reciprocal basis are obtained from 
BP = g«%: 



£ r = Vr = n0) , 
£ s = Vs = — t , 



The gradient and divergence read 



V = & a d n 



(15) 



V-A=-^d ff (V^%) ■ 



(16) 
(17) 



Note that on the interface, the covariant component A r is 
equal to the normal component A„, but that A s is related 
to the tangential component A, by A s = yfgA,, because 
& s is not normalized to one. 

In the following, we will denote inner quantities by 
the uppercase letter corresponding to the lowercase letter 
denoting the outer quantity, whenever this is meaningful. 

Since the interface will move in general and the coor- 
dinates r and s are defined with respect to the interface, 
there is also a transformation rule for the time derivative: 



d,f(x,z, t) = d,F(r, s, t) - vVF(r, s, t) . 



(18) 



where \(s, t) is the interface velocity. Equation ( 18 1 ex 



From (11), we first obtain the metric coefficients g a p — 



hibits that the time derivative in the comoving frame is a 
material derivative. In order to formulate the matching 
conditions concisely, we will occasionally also write the 
outer fields as functions of the variables r and s (without 
changing their naming letter, thus in this case adhering to 
the physicists' convention of using a letter for a physical 
quantity rather than a mathematical function) . 
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B. Inner equations 

To render the scales of the different terms more visible, 



we rewrite Eqs. (16) and (17 1 in the form 



V = ^ n,9 P + ~F t<9s ' 



V ■ A 

y/g — I + WpK 



-t-u 



d p y[g~A r + d s — A s 



(19) 

(20) 
(21) 



Assuming, without loss of generality, that the tangen- 
tial velocity of the interface vanishes, Eq. ^ takes the 
following form 



d,® - i v„5 p <D =V-MV^ 6fi(V 2 <S>, *) , 
5/i(V 2 0,<D)=-^-5 P VI5 p (D 



(22) 



- W £ 



d s O + 2f'(0). 



with 
V ■ MV 



(23) 



Hence, the leading term of the inner equation ( 22 I with 
the differential operator given by (23) is of order W~ 4 . 



C. Expansions, matched asymptotic analysis 

To solve the outer and inner equations successively, 
we expand the phase field in both the outer and inner 
domains in powers of W 

(p(x,z,t) = ^°\x,z,t) + Wcf> m (x,z,t) 

+ W 2 <p {2) (x,z,t)... (24) 

and 

®(r,s,t) = & 0) (r,s,t) + W® m (r,s,i) 

+ W 2 <f> (2 \r,s,t)... . (25) 

We now proceed solving the outer and inner equations 
order by order. 



1. Leading order 

The leading-order outer equation for <p is 

V • MV/'(0 (O) ) = , (26) 

which is to be supplemented with the boundary condi- 
tions 0*°' = 1 and <fP^ — at infinity in the regions where 
the system is solid and non-solid, respectively. If we re- 
gard ( 26 1 as a partial differential equation for the function 



/'(0 <O) ) (rather than for </> (0) itself), this boundary condi- 
tion translates into /'(0 (O) ) — > as [r[ — » oo, which may be 
seen immediately from the explicit form of /'(<A), given 
in App. [B] The new boundary condition is valid every- 
where at infinity except possibly in a region with size of 
order W. For general M{(p, WV<p), the partial differential 
equation (26) is of course nonlinear. Nevertheless, it can 
be shown to have the unique solution /'(0 (O) ) = 0, if M is 
positive everywhere, except possibly on a set of measure 
zero. 



To see this, multiply Eq. (26 1 by /'(0 (O) ), integrate over 



which lim 



all of space and use Gauss's theorem: 

= - J dVM [V/'(0 (O) )] 2 + O(W) , (27) 

where the 0(W) stands for the surface integral at in- 
finity. (We will often use the nomenclature for three- 
dimensional systems, not distinguishing between surface 
integrals and boundary contour integrals; also we em- 
ploy dV for the volume element of both two- and three- 
dimensional space.) If M is positive almost everywhere, 
we immediately get /'(0 (O) ) — const., and the boundary 
conditions require the constant to be zero. This con- 
clusion remains of course unchanged, if M becomes zero 
only when (O) is zero or one - a standard choice [2H] is 

M(0)oc0 2 (l -(f,) 2 . 

Hence, the unique solution to the leading-order outer 
problem is, if we now consider it an equation for 0® 
again, (O) = 1 in Q and <pf® = in where f2 T are 
those regions of space, separated by the interface (s), in 
= 1 and 0, respectively. The solution 
= |, still possible for the equation interpreted as an 
equation for f'((p^), is excluded by the boundary condi- 
tions for (f>®\ (This argument presupposes that we have 
no domains that are not connected with infinity. For the 
interior of a closed interface, the solution (O) = 1/2 would 
have to be excluded by a stability argument or by making 
reference to initial conditions.) 

It is then seen by inspection that the outer equation is 
indeed solved to all orders by the solution under discus- 
sion. With the usual construction of coupling terms to, 
say, mechanical or electrical degrees of freedom [TTJ 127] , 
this remains true for the full phase-field model including 
the coupling, as these terms are typically made to vanish 
in the bulk. So our statement, which has some impor- 
tance as we shall see, is valid beyond the oversimplified 
"free" model considered here for the purpose of demon- 
stration. Any deviation of the outer solution from <p^ 
must be transcendentally small. 

Therefore, we have <jP^ = 0, <p <r> = 0, which provides us 
with partial boundary conditions for the inner solutions 
(£(1)^ <j)(2)^ an( j g0 on ( gee App. [X]). Moreover, only the 
inner problem needs to be considered beyond the leading 
order. 

Because g — 1 + 0(W), the leading-order inner problem 
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becomes [see Eqs. (p2| and (23)] 



3 p M(0> (0) )d p [6> pp O (0) - 2/'(cD (0) )] = , (28) 



where d pp = 3 p . This can be integrated once to yield 



2/'(0> (0) )] 



M(0>(°)) 



(29) 



where ci(s) is a function of integration. It is here that 
we have to follow different lines of arguments, depending 
on whether M approaches zero for p — > +oo, which is the 
case for the mobilities assumed in [25] [26] [28] , or whether 
it is just a bounded (and possibly constant) function of 
<p. In the first case, we may immediately conclude c\ - 0, 
because the right hand side of (29) must not diverge. In 



the second case, we obtain the same result by integrating 
(29 1 first and invoking the boundary conditions: 



«9 pp <D (0) - 2/'(O <0) ) = Cl (s) f P i r dp + c 2 (s) . 

Jo M 



(30) 



Since M is bounded from above and positive, the inte- 
gral will be larger in magnitude than 1 / (sup p M)dp - 
pi sup p M, so the two factors multiplying c\ and c 2 are 
linearly independent. The left hand side approaches zero 
for p — > +oo [the argument will be made more rigorous 
below in the discussion of <5 (1) ], so both c t and c 2 must be 
equal to zero. To argue that c 2 is zero in the case where 
M — > for p — > +oo, we can proceed the same way, except 
that we have already gotten rid of the term containing 
C\, so the right hand side of (30 1 is c 2 only. 



To summarize, the leading-order inner equation results 



3 PP <D (0) - 2/'(O (0) ) = , 



(31) 



and this provides us with the solution <£ (0) = | (1 -tanhp) 
as is shown in App. [33] 



to obtain the same result (where for once we have denoted 
an outer quantity by a subscript "out") . 

Integrating once more with respect to p and writing 
out 5ff x \ we have 



6ft (l) = -3 pp O (1 '-^ p O (0) +2/"(<l) <0) )OW 
= d 2 (s) . 



(35) 



Up to this point, there is agreement between this and 
preceding asymptotic analyses [27], if not in all details of 
the procedure, so at least in the results. 

Let us now try to determine the function of integra- 
tion d 2 (s). A priori, there is no reason to use a procedure 
different from what we have done before. We know the 
limiting behavior for p —* ±oo for two of the four terms on 
the right hand side of ( 35 1 : lim p ^ ±00 <9 p O (0) = [which fol- 



lows from either the matching conditions or by inspection 
of the solution ( B8 1] and lim p ^ ±00 d 2 (s) = d 2 (s) (because 
d 2 is independent of p). Moreover, from the matching 
conditions, we obtain the limit for <S) m 

®m „ p0'(°>(+o) + (I) (+o) = <p m (+0) = o 

(p +oo) . (36) 

The second equality follows from the fact that <^ 0) = or 
(O) = 1, hence its derivative with respect to r vanishes 
on both sides of the interface; the third equality is a 
consequence of the fact that (O) solves the outer equation 
to all orders and hence <^ (1) = 0. In addition, it can be 
shown 29J that if the interaction with mechanical degrees 
of freedom is included, this will only lead to terms that 
also vanish in the limit p — > +oo. 



With three of the four terms in ( 35 1 having a definite 



limit, we may conclude that the fourth must have a limit 
as well and obtain 



lim -d pp 0> (1) = d 2 (s) 



(37) 



2. Next-to-leading order 



The next-to-leading order in Eq. (22 1 is the order W 



Since the differential operator in front of the chemical 
potential is of order W~ 2 and the chemical potential mul- 
tiplied by another factor W~ 2 , we must expand 6fi up to 
order W. Equation (31) already tells us that <5ju (0) 
we obtain 



■ 0, so 



<9 p M(O (0) )<9pS£ (1) = . 



from which we get 



(32) 



(33) 



M(d>(°)) 

As before, we can immediately conclude from this that 
d x = 0, if we assume M(O (0) ) -> for O (0) -> 0, 1. For arbi- 
trary but bounded M, we invoke the matching conditions 
(see App. [A} 



lim tdpdfP = drdfiZu = 



(34) 



But if this limit exists, it cannot be different from 
zero: transforming to £ = 1/p, we see that d pp ^> { ' = 
(fdf) 4> (I) , which impl ies the asymptotic behavior (1) ~ 
-d 2 /2£ 2 — > 0) and hence the divergence of <£ (1) as 
-d 2 p 2 /2, if d 2 + 0. (The same kind of argument can be 
used to show that the left hand side of Eq. ( 30 1 goes to 



zero, even though the matching conditions do not provide 
a direct expression for lim p ^ ±00 <9 pp <t> (0) .) 



The upshot of these detailed considerations is that 



d 2 (s) = 6fi m = 



(38) 



Previous treatments of the problem did not enter into 
these considerations. Instead, one of the following two 
equivalent approaches was chosen. Either, Eq. (351 was 



interpreted as a linear inhomogeneous differential equa- 
tion for ' and Fredholm's alternative invoked. Since 
the appearing linear operator 



£ = 3 pp -2/"(O< >) 



(39) 
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is hermitean, we know (from taking the derivative of 
Eq. § w.r.t. p) that d p O (0) is a solution to the ad- 
joint homogeneous equation. The inhomogeneity of the 
differential equation must be orthogonal to this solution. 
Or else, Fredholm's alternative was not mentioned, the 
equation was simply multiplied by d p O (0) , integrated, and 
it was shown via integration by parts that the terms con- 
taining O' 1 -* disappear. Of course, this is the same thing. 



We then obtain from (35 



k(3 p O (0) ) dp= d p ¥ 0) d 2 (s)dp = -d 2 (s) . 

CO %J — CO 



from which we get, using (BIO I 



d 2 (s) 



1 



(40) 



(41) 



From Eq. (J9|, we gather that an expansion of 6p out in 
powers of W will contain three types of terms, the first 
of which have the form V 2 (f> ( ^ (k = 0, 1 , . . . ) , while the 
second contain factors (ft® (k = 1,2,...), coming from 
an expansion of f'(<p) about (p <0 \ and the third include 
alone. All of these terms vanish, because (p^ = 
for k > and because /'(0 (O) ) = 0. This is simply a 
consequence of the fact that the outer equation is solved 
exactly by <O) = and (O) = 1. The "chemical potential" 
appearing in the phase-field equations needs to be related 
to the true, i.e., sharp-interface chemical potential only 
at the interface. In the outer domain, it is zero. We 
can then conclude from ( |44| that e\(s) = (of course 
e 2 (s) = 0, too, but we shall not make use of that result). 

Given that SpP^ is independent of p, the inner equation 
at order W~ l takes the form 



Both Eqs. (38 1 and (41 1 were derived by valid meth- 
ods, therefore they should both hold true. Nevertheless, 
as we shall see shortly, Eq. (38) is a quite undesirable 
result. This may be the deeper reason why it was so 



far overlooked and only the analog of Eq. ( 41 I derived 



When Eq. (38) is inserted in (41 ), it leads to zero curva- 



ture at lowest order. In a phase-field model for the ATG 
instability the same kind of reasoning imposes a re- 
lationship between the elastic state of the material and 
the curvature. In models, where the interaction term 
is quadratic in W |27j . it again imposes the restriction 
k — 0(W). To summarize, in all cases we obtain a restric- 
tion on the curvature at lowest order, which means that 
the phase-field model will not be asymptotic as long as 
the deviation from this imposed value is not small. 



3. Higher orders 



To see that the model would indeed work if we did not 



have the restriction ( 38 ) , let us consider the equations at 



the next two orders, ignoring for the time being the result 
d% = 0. Since both 6p (0) {= 0) and 6p (l ^ are independent 



of p, the first term of the operator (23 1 does not produce 
any contribution from these terms in (22 1, and the order 
W~ 2 equation reads 



d p Md p 6p (2) + d s Md s 6p (u> = 



(42) 



where we can immediately drop the second term, because 
of = 0. After two integrations this becomes 



Sjj 



r,(2) 



rP ! 

ei(s) -rrdp + e^s) 
Jo M 



(43) 



If M — > for p — > +oo, we immediately find e\(s) — 0. 
In the general case, we use the matching conditions [see 



(A6I] 



6fl 



^p 1 d rr dfif^\ r=± o + pd r 6p^\ r=±0 

~(2) 

+ ^outlf=±0 ■ 



(44) 



- v„(9 p <D (0) = d p Md p 6fP ) + d s Md s dp 



(45) 



After integration over p (y n does not depend on p) we 
find 



Xco 
M(<5 <0) ) dp + M(<D (0) ) dpdfPA 
CO 



.(46) 



Here we can drop the second term on the right hand side, 
if lim p ^ ±M M(<5 (0) ) = 0. Formally setting M(O (0) ) dp = 



3M and using ( |41[ ), we arrive at 

V„ = MdrcK . 



(47) 



Hence, (47 1 reproduces the sharp-interface limit d3l, with 
the relationship between M s and M defined in dTT! 

Finally, M would be infinite for positive functions 
M(O (0) ) that do not reduce to zero for p — > +oo; there- 
fore, in the end we would indeed have to require M(O (0) ) 
to decay far from the interface, if 6fr ' were different from 
zero. In reality, we do not just have (47 1, the equation we 
want, but in addition Eq. (38 1, requiring 6p (l) — and, 
consequently 



v„ = . 



(48) 



At first sight, Eqs. (47 1 and (48) may look like contra- 



dicting each other, as we can prepare an initial state with 
arbitrary curvature of the interface, and hence the veloc- 



ity should be different from zero according to (|47|) but 
equal to zero according to 



However, in preparing 
an arbitrary initial state, we have no certainty that the 
system will already follow its (lowest-order) asymptotic 
dynamics. A similar phenomenon happens in all phase- 
field models when a simulation is started with an initial 
interface perturbed by white noise. Since the asymp- 
totics of the phase-field equations require curvatures to 
be smaller than 1/W, the initial stage of the dynamics 
where larger curvatures are present, will not be governed 
by these asymptotics. But in that case the asymptotic 
behavior is sufficiently robust to keep the initial stage 
short. 
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Since ( 48 1 is a lowest-order result for the interface ve- 
locity, we can conceive of two different scenarios. Ei- 
ther the next-order result for v„ is nonzero. Then the 
model might asymptotically reproduce the desired sharp- 
interface limit at the next order, but its utility would be 
restricted as it would be quantitative only for k — O(W) 
(i.e., kW <s W/L, where L is a typical system length 
scale such as, for example, the equilibrium diameter 
of a crystal). Its validity would be restricted to near- 
equilibrium situations. Or else v„ is zero at all orders of 
the asymptotic expansion, hence transcendentally small 
in the asymptotic limit. Again, this could describe a 
near-equilibrium situation at best. Moreover, such a 
model would violate the spirit of phase-field models in 
general, in which we seek to have an analytic statement 
about the sharp-interface limit with as short an expan- 
sion as possible. The necessity to perform asymptotics 
beyond all orders should not arise in problems where we 
have such a high degree of freedom in constructing the 
model equations. 

As to the general behavior of the SM model, it is quite 
tempting to speculate that when conditions are such that 
(48 1 does not hold yet, the phase-field model discussed 



here will satisfy all the other less restrictive conditions 



already, including (47 1. Then the model would be appli- 



cable during the period where the influence of condition 
(38 1 leading to (48 1 is still small. However, it should 



be clear that without a theoretical estimate of the er- 
ror in this not fully asymptotic state, the model can 
hardly be considered quantitative. Condition (381) should 



be expected to have a stabilizing influence on decaying 
modes of the interface, accelerating their relaxation to- 
wards equilibrium. Its effect on growing, i.e. unstable 
modes is difficult to assess. 

It is instructive to note why the nonconservative model 
obtained when ^ is replaced with (JsJ) does not suffer 
from a similar difficulty. In that model, the velocity is 
already determined at the next-to leading order. Instead 
of ( 35 ) , we get 



Vn d p ^ = K\d pp & l) + Kd p & Q) - 2/"(0<°W 1 > 



(49) 



Again we may conclude that all the terms on the right 
hand side go to zero as p is sent to +oo. However, this 
does not lead to any constraints, since the left hand side 
is p dependent now and goes to zero as well, satisfying the 
limit automatically, whereas in the surface-diffusion case, 
it was a function of s only (d.2) that could be concluded 
to be equal to zero. So consideration of the limit does 
not produce anything new here, and the only procedure 
available to extract information on v„ is to use Fredholm's 
alternative which gives the correct sharp-interface limit. 

In the case of the nonconservative model, the intro- 
duced chemical potential functional is zero in the bulk 
just as in the conservative case, but there are no restric- 
tions on its variation near the interface, where it acquires 
a form tending to a 5 function in the sharp-interface limit. 



In the conservative model, this is excluded by restrictions 
on the derivative of the chemical potential with respect 
to p, meaning that the latter must be smooth across the 
interface. Since it is zero off the interface, it is zero on it 
as well. Due to this reason, the phase-field model strictly 
speaking applies only to the equilibrium limit. Far-from 
equilibrium dynamics is not likely to be captured faith- 
fully. 

Out of the phase-field models for surface diffusion con- 
sidered in the literature, the only one that is (apart from 
our own work [29 ) not subject to the criticism offered 
here seems to be the one given by Ratz, Ribalta, and 
Voigt [23] . Let us briefly discuss the asymptotics of this 
model that we will henceforth denote as the RRV model. 
In their simplest form, i.e., for isotropic surface tension 
and vanishing kinetic coefficient, the model equations 
read 



dep 
~3t 



= Vj 



j =MB(0)V— <5//(V 2 <M), 
^ = ^)(- W2V20 + 2/ '^)' 



(50) 



with mobility function = 120 2 (1 - </>) 2 , double-well 
potential f(<p) = 2 (1 - 0) 2 , and the so-called stabiliz- 
ing function g{<p) = 1O0 2 (1 - (p) 2 . Here, we have rescaled 
the equations from [21] so as to obtain the same zeroth- 
order interface profile as in the SM model (with the orig- 
inal equations, the interface would have one third of the 
width ouf our profile). The leading-order inner problem 
becomes (311 again. At next-to leading order, we obtain 
<5/z (1) = da(s) (as before), but now the chemical potential 
function is defined differently - it has a prefactor that 
diverges in the bulk 

s» m = ^(-^a»w-^°) + 2/''(^)a>™). 

(51) 

(The first-order term due to variation of the denominator 
vanishes, as it contains the differential expression from 
the left hand side of (311 as a factor.) The numerator 
of the right hand side of pTj ) goes to zero as p — > ±00 
but so does the denominator g(<£> (0) ), which renders 5p m 
indefinite, thus introducing the degree of freedom neces- 
sary for a nonzero value d2(s). Multiplying the equation 
by g(cD <0) )<9p <E> (0) and integrating with respect to p from 



-00 to 00, we arrive at 



f 

•J — 1 



d 2 I 8 (® (0> )d p & u> dp 



f.(0). 



00 



dp , (52) 



where use has been made of the fact that <5 (0) is a left null 
eigenvector of the linear operator [inside the parentheses 



in J5T |] acting on (1) , to get rid of the O a) terms. The 
integrals in (52) are evaluated in App.jB] they are equal 
to -1/3 and 1/3, respectively. Hence, d 2 — k. 
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The steps for the following two orders of W follow pre- 



cisely the scheme leading from (42 1 to ( 47 1 , which then 
yields v„ = M*d ss K, where M* = J^MB(O (0) ) dp = M (for 
the integral see App . |B| , hence we obtain the desired 
sharp-interface limit (|3|. 

Note that with this model, it is essential that the mo- 
bility function goes to zero off the interface. For the 
chemical potential dp varies in the bulk (it behaves as 
d2(s) near the interface), hence diffusion there must be 
suppressed by a vanishing mobility. 



IV. TENSORIAL MOBILITY 

While the RRV model avoids the mistake of imposing 
(38 1, it does so by a purely mathematical device, the 

It is then 



introduction of the stabilizing function g((f>). 
natural to ask whether an accurate model may not be 
derived on the basis of mainly physical considerations. 

That the phase-field model given by Eq. (J9j) does not 
quite yield the correct asymptotics may be traced back to 
the fact that the differential operator V ■ MV, prepended 
to the chemical potential, does not reduce to the surface 
Laplacian A s in the asymptotic limit. In fact, the second 



term on the right hand side of Eq. (23 \ is, up to a factor, 
the Laplace-Beltrami operator on the surface (for p — 0) , 
but the first term, containing derivatives with respect 
to p is orders of magnitude larger, being preceded by a 
factor of l/W 2 . As a consequence, the asymptotics must 
be secured by the full solution of the equation rather 
than by both the operator and the chemical potential 
converging to the desired sharp-interface limits. 

Realizing this property of the model, it seems natural 
to modify the differential operator via introduction of an 
essentially tensorial mobility. Let us denote by 



n = - 



V0 



(53) 



the normal on the surface <p = const, (for (p — 1/2, we have 
n = n) , then we expect the operator V • PV with 



1 



n : n 



(54) 



(cartesian components: Pjj — Sij — hfij) to reduce to the 
surface Laplacian asymptotically. A colon is used to de- 
signate a dyadic product, so P is a projection operator 
projecting onto the tangential plane of a level set of (p. 
Introducing the shorthand <f>, a — d a <p, we have V<p — & a (f>, a 
and 

V(l-fl:ft)V = 

^ yfgg" [d, - ty) ■ (55) 

The third expression is obtained from the second apply- 
ing the divergence operator (17) and renaming a — * p, 



P — > a, p — > p in the three pairs of "mute" indices. 



To expand this operator in powers of W in the inner 
domain, we introduce an abbreviation for a normalized 
gradient of <£, being of order W°: 



(va>) 2 s ^(vo>) 2 , 

(V<D) 2 = W 2 ®, a g a ^,B = 0> 2 + 



w 2 



(56) 



-<D- 



Inserting this into ( 55 1 and carrying the expansion to 
formal order W , we find first that the order W~ 2 terms 
(containing two derivatives with respect to p) cancel each 
other. The remainder reads 



V-(l-ft:ft)V = 



1 



8 n 



1 



(® 2 d p -®, p <!>, s d s ) 



<» v o„. 

(V<D) 2 



d p \ + 0(W), 



(57) 



and this expression still contains derivatives with respect 
to p. However, if the leading-order solution <D (0) depends 
on p only, as it did in the last section, then all the deriva- 
tives of <I> with respect to s are 0(W) at least, and since 

= 1 + 0(W), Eq. ((57) reduces to V • (1 - n : n)V = 
<9 2 + 0(W), i.e., at leading order the operator indeed be- 
comes the Laplace-Beltrami operator on the surface. 
This then suggests to replace the phase-field equation 

with 

d -l = MV • (1 - n : n) V T^Sp(V 2 cf>, 0) , (58) 
at W £ 

where dp is unchanged from ((£]) but M is a constant mo- 
bility now. 

In this model, the equation for the velocity would ap- 
pear at the next-to leading order already and take the 
form 



v„d p o (0) = Md s 



5 PP <D (1) + /«9p<D (()) - 2/"((D (0, ) <£ 



(59) 



Because the operators X [defined in Eq. (39 1] and d ss com- 
mute, Fredholm's alternative is applicable the same way 
as in the nonconservative case. This eliminates O fl) from 
the equation and produces the correct sharp-interface 
limit . 



In spite of this enjoyable state of affairs, model (58 1 
fails much more miserably than ([9j. The reason is that 
the zeroth-order solution is not unique. In fact, the 
leading-order outer equation 



V-(l-n< o >:n (O) )v/'(0<°>) = O 



(60) 



is solved by any (differentiable) function (O) satisfying 
the boundary conditions: obviously we have V/'(0 (O) ) = 
/"(<^ (0) )V^ (0) , whence 



A (0) . fl (O) V /'(0W) = -ftW/"(0W)|V0' 



x(0) 



(0), 



a(«)i 



= f(^W (0 » = V/'(f), (61) 



(OK 
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which implies (1 - n (0) : n (O) )V/'(0 (O) ) = for all func- 
tions /' of <O) . Intuitively, this behavior can be easily 
understood for a planar interface. Then the equation 



of motion (58) strictly contains only derivatives of the 
phase field parallel to the interface, and the profile in 
the perpendicular direction therefore remains completely 
undetermined. 

It can be said that this model fails for reasons comple- 
mentary to those of the scalar model. Whereas we had 
one equation too many in that case, adding a constraint 
to the desired sharp-interface dynamics, now we have one 
equation too few, as there is nothing in the model fixing 
<p a)) . If we had the right (O) , the tensorial model would 
work perfectly. 



V. MODIFIED TENSORIAL MOBILITY 
MODELS WITH CORRECT ASYMPTOTICS 



In order to obtain a model not plagued by either of the 
disadvantages of the two cases discussed, it appears that 
it is useful to combine ideas from both. While it is cer- 
tainly desirable to have a differential operator that itself 
approaches the surface Laplacian, it should do so only for 
phase field functions that have the correct leading-order 
profile. 



These considerations lead us to make the ansatz 
MVj 



dcf> 
~dt 



1 



j = e m v— 

5fi = -W 2 V 2 + 2/'(0) 



(64) 



and leave the precise choice of the value of m for later 
it will be suggested by the asymptotic analysis. 
The corresponding inner equations are 



w w z 



<5ju(V 2 1>,<t>) 



(65) 



-W z 



with 



W 2 (V<D) 2 
4/(0) 



(VO) 



(66) 



Leading order 



A. Locally conservative model 



One way to achieve this goal is to modify 1 - n : n into 

, V0 : V0 W 2 CVd>) 2 , , 

Q=l-W 2 * rt , = 1 - - - ' n : n ■ (62) 



4/(0) 



4/(0) 



If we replace the projection operator in Eq. ( |58| by 
Q, then the outer equation at leading order will have 
the same differential operator as the scalar model with 
constant M. 

On the other hand, in the inner domain, we have, 
provided cD (0) solves the differential equation ( B6 1 , 
W 2 (V0>) 2 = Ol+OiW 2 ) = 4/(<D) + 0(W) [this follows from 
Eqs. @, flB7| , and ([Bl])] , whence Q*l-n:n. 

This approximation is accurate up to 0(W) only, which 
is not sufficient, because the order W correction would en- 
ter as a bothersome additive term in the next-to leading 
order inner equation. 

A better inner approximation to 1 — n : n than just Q 
is provided by a minor modification. Obviously, we have 
Q- 1 - n : n + <3(W) n : n in the inner region. Taking this 
to some integer power m we get, because 1 — n : n and 
n : n are orthogonal projectors: 



Q m = 1 - n : n + 0(W") n : n . 



(63) 



In the outer equations, Q becomes the identity opera- 
tor to leading order, i.e., (2 <O) (0 (O) ) = 1, and at the lowest 
order in W, we have 



V 2 /'(0 (O) ) = 



(67) 



a Laplace equation that we know to be uniquely solvable 
for /'(0 <O) ) with Dirichlet boundary conditions at infin- 
ity. This boundary condition is even homogeneous (ex- 
cept possibly in a part of the boundary region at infinity 
having a size of order W) , leading to the unique solution 
/'(0 <o) ) = 0. This leaves us with the three possibilities 
(O) =0, \, 1, of which (O) = or (O) = 1 are realized, 
according to the particular boundary condition on <p®\ 

Again, = and 0=1 are solutions to the outer 
problem at all orders of W. Admittedly, the operator 
Q becomes indefinite at order W 2 for = and = 
1 [because of the denominator /(0)], but this does not 
matter, since the expression for 6jl alone is zero already 
at = and 0=1. 

The leading-order inner equation reads [g = 1 + <9(W)] 



1 - 



(0) y 



4/(0<°)) 



d p( d PP® (0) ~ 2/'(1> (0) )) = . 



(68) 



Clearly, this is solved by O (0) = ^ (1 -tanhp), which makes 
both the expression in brackets and in large parentheses 
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vanish. If we require m to be even, this solution is more- 
over unique (up to translations, which are eliminated by 
the requirement that the interface be at p = 0). For as 
soon as we assume (<D ,f> f t 4/(<5 (0) ), the mth power of 
the bracket expression will be positive, allowing us to use 
similar arguments as in Sec. IIIC between Eqs. (28) and 
Q to prove that d pp & ()) - 2/'(<D (0) ) = 0, and hence the 



bracket expression must be zero, contrary to our assump- 
tion. Thus we do get a definite solution for (D (0) from the 
inner equation, which moreover shows that at leading or- 
der of the inner expansion the second-order p derivatives 
of the operator V • Q"'V cancel each other. 



2. Next-to-leading order 

To simplify computations at the next order, we first 
expand V ■ Q'"V up to formal order W°. This produces 



V ■ g'"V = 



w 2 cv<&y 



4/(0) 



i 



1 <D, 



~F d p—p^ d * F "s—p-z—Op 



1 



1 <E>2 



(69) 



Given that (D (0) is a function of p only, we realize that the 
third and fourth terms on the right hand side are O(W), 
containing derivatives with respect to s of <£, the fifth is 
even 0(W 2 ), so these terms may be dropped immediately 
in an expansion up to 0(1). The second term on the 
right hand side owes its existence to the fact that Q is 
not exactly the projection operator on n [note that no 
such term is present in Eq. (57l] and it has a prefactor 
of 1 /W 2 due to the double derivative in p. This term 
which is desirable at leading order, because without it 
we would not have a determinate zeroth order solution 
O^ \ is somewhat disturbing at the next order. Since the 
order of this term is 0(W m ~ 2 ), we can make it small by 
choosing m > 3, i.e., restricting ourselves to even m for 
the reasons discussed before, we set m — 4. Then the only 
remaining term on the right hand side of Eq. ( 69 1 up to 
order W° is the first term, which is the desired surface 
Laplacian. 

Using this result, we can write the next-to- leading 
(nontrivial) order inner equation 



= Md ss 6p w 
6p w = 



(70) 



again with X as given in Eq. ( 39 I 



Note that we actually seem to have skipped orders 
here. The leading-order inner equation is formally 
0(W~ 4 ), but once the zeroth-order inner solution is fixed, 
the differential operator V • Q'"V is, according to (69 1, of 



order y^ ma x((),m-2) on i Vj so t ne orc j e r W~ 3 vanishes identi- 
cally. The order W~ 2 is satisfied automatically, because 
the zeroth-order chemical potential is zero; the next non- 
trivial order is W~ l . Alternatively, one may say that the 
effective leading order has become W~ 2 . 

The total linear operator in front of <E> (1) becomes 
—d ss £.. It is hermitian, because its hermitian factors com- 
mute. Hence, <9 P <I> (0) is a left null cigcnfunction. Multiply- 
ing ( 70 I from the left by it , integrating with respect to p 
from -oo to oo, we obtain Eq. (47 1. This proves that the 



considered phase-field model based on a modified ten- 
sorial mobility has the correct asymptotic behavior for 
small W, neither overconstraining the system by adding, 
nor leaving it indeterminate by losing equations. 

Clearly, Eq. (64) establishes a local conservation law 
for (p, i.e., the rate of change of the integral of <p over 
some control volume is given by the integral of the cur- 
rent j associated with cp over the surface of the volume, 
and this holds for arbitrarily small volumes. <p is the den- 
sity of a conserved quantity. In particular, for a system 
with boundaries through which there is no flux, the vol- 
ume integral of <p will be conserved. Therefore, we will 
denote the model discussed in this section as the locally 
conservative tensorial or LCT model. 



B. Globally conservative model 

If one is willing to give up the conservative nature 
of the phase-field equations themselves and requires the 
conservation law only in the asymptotic limit, an even 
simpler construction is feasible. 

Consider the model 



^ = MV ■ PV^Sp+N 
dt W 2 H 



(n-V) 2 0-— /'(0) 



1 



n : n 



(71) 



6fi = -W 2 V 2 4> + 2/'(0) 



with both M and positive constants. With N = 0, this 
reduces to the tensorial model of Sec. |IV| With jV > 0, it 
can be shown along similar lines as in Sec. Ill C 1 that 



the leading-order outer solution for dp is unique leading 
to the solutions cf> — and <f> — I, depending on the bound- 
ary conditions at infinity. Moreover, the leading-order in- 
ner solution with boundary conditions linip^-oo <t> (0) = 0, 
lim p ^ M <J> (0) = 1 can be shown to be unique up to trans- 
lations along p and is given by (B8) after requiring p — 
to correspond to the value \ of the phase field. 

The role of the term is only to fix the profile of the 
phase field at leading order, otherwise it is constructed 
so as to not affect the normal velocity of the interface. 
Once is set, the next-to- leading order inner equation 
reads: 



(Md ss -N)£® a> = v„«9 p (D ( 



f,(0) 



Md ss Kd p (t> (0> 



(72) 

and since Md ss — N commutes with X, we obtain the de- 
sired sharp-interface limit again. Our numerical investi- 
gations indeed show that the results depend only weakly 
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on the choice of the parameter N, even for moderate sep- 
aration of the length scales. We find N < 2.5M/W 2 to 
already give satisfactory results - there are only small 
differences to results obtained when N is two orders of 
magnitude smaller, i.e., for N = 1.25 x 1(T 2 M/W 2 . 



While the model (71 1 is asymptotically conservative, 



it is desirable to have exact global conservation of the 
phase-field, because this will render long-time simula- 
tions more reliable. As the model stands, one might be 
obliged to choose the interface width smaller as the simu- 
lation time becomes larger, which is certainly something 
one would wish to avoid. Therefore, even though the vi- 
olation of phase-field conservation is small and the model 
would already be useful in its present form, let us look for 
an improvement restoring global conservation. By this 
we mean that <j) need not be the density of a conserved 
quantity, hence its time derivative need not be the diver- 
gence of a current, but for no-flux boundary conditions, 
the total volume integral of (p should remain conserved. 
This can of course be achieved via the introduction of a 
Lagrange parameter: 



at 



= MV 



W 2 H 



(n-V) 2 -—/'(<£) -A(r,0. 

(73) 

Here, we have allowed A to depend on r which gives 
useful additional freedom for improvement of the model 
as we shall see immediately. If A were restricted to being 
a simple number, it would have to to have the value 



A = 



vj v 

X 



dV 



d<f> 



Id 



at 



dV 



(ft ■ V) 2 



W2- 



(74) 



where d<f> \d/dt is the time derivative of the phase field ac- 



cording to ( 71 1 and V is the volume (or area, in 2D) of the 



system. Since the first term of the right hand side of ( 73 ) 



is conservative anyway, it drops out of the calculation of 
A, if no fluxes through the system boundary are present. 
A drawback of the formulation ( 74 ) is that it would lead 



to a modification of the phase field in the bulk from the 
equilibrium values and 1, as soon as the Lagrange mul- 
tiplier became nonzero. This can be avoided by taking 
advantage of the liberty to make A vary in space (i.e., 
we consider a whole set of Lagrange multipliers, not just 
one). If we take A of the form 



A(r,r) 



|V0| 



j v dvm 



X 



dV 



(ft ■ V) 2 



W 2 



(75) 



the global conservation law is restored without any mod- 
ification of the bulk solutions. We will call the model 



described by Eqs. (73) and (75 1 the globally conservative 



tensorial or GOT model. 

Of course, we have to verify that the introduction of 
the Lagrange parameter does not destroy the asymptotic 
validity of the model. Clearly, the parameter disappears 



from the leading order of the equation; but the inter- 
face velocity is determined at next-to leading order, and 
in general one would expect A(r, t) to contribute to the 
equation at that order. This turns out not to be the 
case and is due to the judicious choice of the form of the 
parameter, as we shall see now. 

The next-to leading order inner equation can be writ- 
ten 



(Md ss - AO + N 



J dVdp&V 



X 



dV£® 



Md ssK d p ¥ {)) 



(76) 

and we are in the awkward situation that the linear oper- 
ator acting on C> (1) is not self-adjoint, so the application 
of Fredholm's alternative becomes nontrivial. However, 
the equation contains several terms oc d p <f> <0 \ which sug- 
gests to have £, act on it, leading first to the much simpler 
equation 

£(N -Md ss )£<& m = 0, (77) 

because -C3 p O (ff) = 0. But here the operator acting on <t> (1) 
is semipositive, the operator sandwiched by the two Xs 
strictly positive. Hence, we may conclude that £<£> w - 0. 
But then the left-hand side of ( 76 ) is zero, meaning that 



the linear equation is in fact homogeneous and the right- 
hand side has to vanish, too. This implies 



v„ = M8 ss k , 



(78) 



the sought-for asymptotic result for the interface velocity 
It also implies that the Lagrange multiplier is 0(1), in- 
stead of 0{W~ X ), i.e., it is by a factor of the order of (Wk) 2 
smaller than the leading-order terms of the equation. 

This supports what we can point out on the basis of nu- 
merical studies: for reasonable separation of length scales 
as they appear in typical simulations, the influence of the 
Lagrange parameter is negligibly small, at least for not 
too long time scales. 



VI. ANALYTIC LINEAR STABILITY 
ANALYSIS 

A linear stability analysis of a stationary planar front 
may be useful in trying to differentiate between the mod- 
els. Clearly, one should require that the spectrum ob- 
tained for the phase-field model reduces to that of the 
sharp-interface model [e.g., Eqs. ^ or ([7])] in the limit 
of small interface width. 

For simplicity, let us first consider the nonconservative 
model described by A planar front solution is given 

by 

= o (z) = <D (Z) = ^(1-tanhZ) , Z = i- (79) 
2 W 

Adding a small perturbation 5<p, i.e., setting <p = <po + 
6<j)(x,z,t), we obtain the linearized equation 

dS<p . 2 



dt 



W 2 



f"(<h)W- 



(80) 
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Using the ansatz 



%z,t) = ^{Z)e 



ikx+iot 



we obtain the eigenvalue problem 



w¥=^ (dzz - W 2 k 2 - 2/"(O )) V . 



w 2 



(81) 



(82) 



While this is a linear problem, it is one that does not have 
constant coefficients, so an exact solution is not readily 
available. Instead, we must rely on asymptotic analysis 
again to make progress. However, as it turns out by 
reinserting the found eigenfunctions and eigenvalues into 
(82 1, the expansion provides exact results in this case. 



Details of the calculation are given in App. [C] We find 
two branches of the spectrum with 



AM , 



(83) 



-Mk 1 



For the first branch, the eigenfunction does not vanish 
as Z — > +oo, for the second it does and moreover is pro- 
portional to Oq(Z). As Wk <K 1, any contribution of the 
perturbation containing the first eigenfunction will de- 
cay fast, leaving a remainder that decays with rate a>b, 
which corresponds to the dispersion relation of the sharp 
interface- limit, Eq. ([6]). Note also that if we assume our 
initial perturbation to describe a slightly perturbed pla- 
nar front, i.e., 



= i{l-tanh[Z-«5<r(x)]}, 



then we have 



d® 
' dZ 



si, 



(84) 



(85) 



hence the perturbation is oc <S)' (Z), so the relevance of the 
second eigenvalue is obvious. 

With these preliminaries, the linear stability analysis of 
the LCT model becomes more or less trivial. Linearizing 
the equation of motion (64) about <f>o, we obtain 



ddep 
~dt 



where 



= MVe»-V 2 +p/"((l.o))#, 



Qo = Q(®o) = 1 - e. : e z , 



(86) 



(87) 



as the level set normals are in the z direction and 
W 2 (V(D ) 2 /4/(<D ) = 1. But then Q 4 = Q , and Vfi*V = 
d 2 . Hence, 

854> 



and after inserting ansatz (81 1, one gets 

aW = ^ [d zz - W 2 k 2 - 2/"(€>o)) ¥ . 



(89) 



But this is the same eigenvalue problem as ( |82| with M 
replaced by Mk 2 . So we obtain a two-branch dispersion 
relation again, this time with eigenvalues 



\Mk 2 
W 2 
-Mk\ 



Mk" 



(90) 



which in the light of the preceding discussion is a rea- 
sonable result (the eigenvalue with large absolute value 
is negative and influences the dynamics for a short time 
only). 

Analysing the GCT model is only slightly more in- 
volved, at least in the form without the Lagrangian mul- 
tiplier. One first derives that, for a perturbation about 
Oo the simplification 



6(iiV) 2 <t>= yfidzzi 



(91) 



holds true. Then the linearization of Eq. (71 1 reads 
86<p 



at 



= Md xx (-V 2 + — /"(O ))<ty 
+ §2 ( 5 zz -/"(<*>())) <50. 



(92) 



After using the ansatz (81 1, one can cast the resulting 



eigenvalue problem into the form (see App. |C| 
(o) - Nk 2 ^ = (dzz ~ W 2 k 2 - 2/"(d» () )) T . (93) 



This again has the same form as (82), now with M re- 
placed by Mk 2 + N and u replaced by w - Nk 2 , leading 
to 



4(Mk 2 +N) , 

(O a = z - Mk 4 



W 2 



(94) 



-Mk 4 



Considering now the form with the Lagrangian multiplier, 
( 73 1 with ( 75 1 , we note that 



SA = 



IVOol 



- f d 3 x^(d zz -2f"m)8cf>, 
'ol Jv W 2 



(95) 



which becomes zero, if we insert the eigenfunction $>' (Z) 
belonging to a>b (see App.[C|. So this eigenvalue, which is 
the relevant one, remains unchanged. We can no longer 
exactly calculate the other eigenvalue nor the correspond- 
ing eigenfunction, but we anticipate that this eigenvalue 
still behaves as 1 /W 2 at leading order and is negative, so 
it does not dominate the asymptotic behavior. 

Unfortunately, an analysis of the SM model along the 
same lines turns out to be much more complicated. This 
may be traced back to the fact that the mobility function 
vanishes in the outer domain, leaving no useful linearized 
outer equation. In fact, the outer equation formally be- 
comes d t 6(p — 0, saying that at linear order perturbations 
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will not decay at infinity. In reality, this means that lin- 
earization is not legitimate, as perturbations will decay 
via nonlinear relaxation the leading-order nonlinearity 
is larger than the (vanishing) linear expression. An anal- 
ysis of the linearized inner equation with the requirement 
that d<$> — > for Z — > +00 does not produce any solutions 
for eigenvalues assumed to diverge more slowly as 1 /W 4 
(or not at all). If the requirement 6<t> — > is given up, 
one may construct perturbation eigenfunctions, but these 
contain polynomially diverging terms at infinity. On the 
other hand, eigenvalues of the form to = to^/W 4 + ... 
could not be investigated by asymptotic analysis, because 
with this assumption the lowest-order perturbative prob- 
lem remains unsolvablc analytically. Note, however, that 
one eigenvalue of the form to = -Mk 4 + 0(W) has to be ex- 
pected for any phase-field model trying to meaningfully 
approximate the sharp-interface dynamics 

It is tempting to speculate that the fact of an unde- 
sirable restriction arising from the asymptotic analysis of 
the SM model and the unfeasibility of asymptotic anal- 
ysis in a linear stability calculation of the model have to 
do with each other. This fits nicely with the observation 
that the linearized RRV model does have a usable outer 
equation. While its mobility function B(cf>) is zero in the 
bulk as well, the presence of the stabilizing function g(tp) 
prevents the reduction to a static result. Writing the 
model in one equation 



dt g((f>) 



W 2 



/'( 



(96) 



we can recast the differential operator in front of the 
parentheses 



1 

6 



9 B(4>) 1 

g(<f>) g(4>) 



„ 7 12 Vtb 



(97) 



and the last term remains regular. The linearized equa- 
tion of motion then becomes 



a*>6 

dt 5 



which for |Z| » 1 turns into 



dd(f> 6 
Ik ~5 



M 



V + w si S nz,9 z 
W 



-v 2 + - 2 f"(%) 



w 2 



(98) 
(99) 



a perfectly sensible outer equation. 

Because the inhomogeneus fourth-order linear differen- 
tial equations resulting from ( 98 1 at successive steps of 



the expansion are difficult to solve (this is in part due to 
the linear operator being nonhermitian) and the impossi- 
bility of getting analytic results for the SM model neces- 
sitates a numerical investigation anyway, we will not pur- 
sue the analytic discussion of the RRV model any further. 
We have checked that <5>' {Z)e lkx+OJt is not an exact solu- 
tion to the linearized equation here, so it appears quite 
likely that the relevant branch of the dispersion relation 
contains W dependent corrections, i.e. to — -Mk 4 + aW 2 . 



VII. NUMERICAL SIMULATIONS 

To determine growth rates numerically, we simulate 
the temporal evolution of a planar front subject to a si- 
nusoidal perturbation. This allows the empirical deter- 
mination of the main branch of the dispersion relations 
of the models considered, i.e., of the largest eigenvalue 
(ojb). Moreover, to assess the behavior of the models 
in a growth situation, we include the coupling to elastic 
fields, i.e., we simulate the Grinfeld instability. This is 
easily achieved by replacing the "free" chemical poten- 
tial 6p with one that includes the correct elastic energy 
contribution. 

In particular, we set 

5fi(V 2 (f>, 4>) = - W 2 V 2 (j) + 2f '{(/>) + ^-h'(tf>)V e i , (100) 



with 



37 



ij z k 



(101) 



where G = £'/[2(l + v)] is the shear modulus or first Lame 
constant, A = £V/[(1 + v)(l - 2v)] the second Lame con- 
stant in the solid phase and G and A are the corresponding 
quantities in the second phase. If the second phase is a 
liquid, G = 0, if it is vacuum, G — A — 0. From now on, we 
will focus on the vacuum case, as it is particularly inter- 
esting for diffusion along a free surface. Equation (100 1 



was originally introduced as a phase-field model for two 
coherent solid phases in contact with each other in [31j . 
Moreover, it has been discussed in [TH] that modeling a 
liquid (or vacuum) as a shear-free solid faithfully cap- 
tures the physics of the system in spite of the property 
of formally coherent strains. This is essentially due to 
the fact that only the divergence of the strain tensor is a 
physically meaningful quantity in the shear-free phase. 

By Ui, we denote the displacement field. In addition, 
equations for the strains My = 5 (diij/dxj + duj/dx,^j have 
to be provided, which are given by the mechanical equi- 
librium conditions for the generalized stress tensor cry 



Z 



dd ~u 
dx. 



= 0, 



o-jj = h(4>) o-jj 



(102) 



Uij via Hooke's law 



The function h{<p) is defined in ( B4 ) and o-jj is related to 
law 

= TT^( M ' 7 + T^2 M »^)- ( 103 ) 



Vi, 



That this model produces the right coupling to elastic- 
ity in the sharp-interface limit has been shown in various 
places, e.g. in [TB] und [23]. A similar modification of 5(i 
leads to the fully coupled RRV model |24j . 

To simulate the Grinfeld instability in a strip geometry, 
a dimensionless driving force F is defined as 



F = 



6 2 (A + 2G) 
4yL 



(104) 
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with 6 being a fixed displacement by which the strip is 
elongated in the direction parallel to the interface. We 
use a square system of length L and uniform grid spacing 
Ax. The interface width, a purely numerical parameter, 
is chosen to be W = 5Ax. For a sinusoidal perturbation 
of a uniaxially strained surface by 6y(x) — Ao sin(kx) we 
use a fixed wavenumber kL — An and a small amplitude 
Ao& = 7r/20. To obtain a good separation of the character- 
istic wavelength of the pattern and the interface width, 
we use kW — 0.16. The imposed uniaxial stress is given in 
the figure caption for each case. After having determined 
the maximum admissible time step, as discussed below, 
we take as simulation time step At — 5 ■ 10~ 4 (Ax) 4 /M 
for the scalar models and At = 5 • 10~ 3 (Ax) 4 /M for the 
tensorial models. The Poisson ratio is chosen to be 
v = A/[2(A + G)] = 1/3. For the GCT model, we use 
N = 1.25M/W 2 , and for the SM model the standard 
choice M(<p) = 36M0 2 (1 - cf>) 2 . To minimize the influ- 
ence of the boundaries, we use helical boundary condi- 
tions in the x direction for the displacement fields, i.e. 
u x (L,y) = u x (0,y) + 6, u y (L,y) = u y (0,y). 
Introducing the Griffith length L c 



Lr, = 



8F(1 -2v) ' 



(105) 



we can write the dispersion relation, i.e., the spectrum of 
the linear stability operator, as follows 



' k 3 

oj = M\ k L 



(106) 



meaning that we have unstable modes at small k (k < 
1/L(j) and stable ones at large k (k > I/Lq). To obtain 
the spectrum numerically, we vary L c in the simulations. 

We use the same algorithm as in [3T] , which even works 
for dynamic elasticity. But since we investigate the be- 
ginning of the Grinfeld instability, the observed interface 
velocities are very small in comparison to the speed of 
sound, and the equations effectively reduce to the static 
elastic case of Eq. ( 102 1. 



While the SM model can be discretized in a relatively 
straightforward manner, some care has to be taken in the 
other models to avoid divisions by zero. This is rather 
harmless in the GCT model, where the problem only 
arises in the computation of the vectors n (|V0| goes to 
zero far from the interface but is positive otherwise, so 
it is sufficient to ensure that the denominator of n docs 
not become smaller than a small positive number). The 
main requirement in the RRV model is that g(<p) should 
not be set exactly equal to zero. 

In the LCT model, more attention has to be paid to the 
situation where f(<p) becomes small, as will be discussed 
below. 

Essentially, we make four types of comparison. First, 
we compare the time evolution of sinusoidal fronts (ini- 
tialized with the correct width of the profile) for a number 
of imposed uniaxial stresses and obtain the linear stabil- 
ity spectrum numerically. Second, we increase the time 



step in the simulation for given mesh size until we reach 
the maximum possible time step providing convergence 
to the correct interface dynamics and then compare the 
achieved values. Next, we initialize a planar profile with 
the wrong interface width and observe relaxation to a 
profile of the correct width. In a realistic simulation, 
slight deviations from the correct profile width may eas- 
ily appear in an initial condition for a curved interface, 
as analytical expressions for const ant- width profiles at 
arbitrary curvature are not readily available (even for an 
initial germ with a shape as simple as an ellipse it is not 
quite trivial to give such an expression) . Any phase-field 
code should be robust against these local variations of 
the profile width and should have it relax to the correct 
value. Finally, we look at the evolution of an elliptical 
inclusion. Since the phase-field parameter is a conserved 
quantity, the ellipse should morph to a circle with the 
same area. 

In Figs.[T]to|3] we show the temporal evolution of a sine 
profile starting with a prescribed amplitude for different 
values of the imposed uniaxial stress. The four models 
are compared directly with the sharp interface prediction 
resulting from Eq. (1061. Fig. [l] exemplifies the stress- free 
case discussed analytically. 
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FIG. 1: Amplitude evolution for a uniaxial stress of F = 0, 
i.e., a Griffith length L a = oo. 

All the situations considered correspond to either weak 
decay or weak growth of the amplitude, as the expected 
exponential behavior still appears linear on the consid- 
ered time scale. 

We note that all the models agree with the predicted 
behavior of the sharp-interface limit to within better than 
one percent for our parameters and time span. While it 
may be observed unambiguously that the SM model dis- 
plays the largest deviation from the desired result, one 
may find it surprising that it reproduces the limit so well 
after all, taking into account that it does not have the 
right asymptotics. Presumably the general idea men- 
tioned after Eq. (48 1 is not too far from the truth: those 
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FIG. 2: Amplitude evolution for a uniaxial stress of F = 2. 



FIG. 4: Full spectrum of the ATG instability. 
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FIG. 3: Amplitude evolution for a uniaxial stress of F = 3.5. 



equations of the asymptotic behavior to which the system 
can adjust locally act as an attractor for the dynamics 
even before the full set of equations, implying more global 
restrictions such as Eq. (48 1, becomes active. It is strik- 



ing that this seems to work even in a growth situation, 
where interface velocities increase on average. 

Figure [4] gives a comparison of the linear stability spec- 
tra, obtained by simulation of the four models, with 



the analytical expression Eq. (106) of the sharp-interface 
model. It is pretty clear that the SM model is far- 
thest off the correct value both below and above the 
fastest-growing wavenumber. The LCT model is good 
for wavenumbers above that of the fastest-growing mode 
but shows stronger deviations than both the RRV and 
GCT model below that mode. The latter two models are 
about equally close to the correct spectrum throughout 



FIG. 5: Allowable maximal time step Af as a function of mesh 
size, where I is an arbitrary length unit. Lines are a guide for 
the eyes only. 

Fig. [5] displays the maximum time step for a given 
grid spacing leading to smooth growth where the results 
for all models agreed perfectly, independent of the time- 
discretization. We observe, for all models, a scaling of 
the maximum admissible time step as At ~ Ax 4 /M, which 
is not too surprising given the fact that the equations 
simulated are fourth order in space and first order in 
time, and we used straightforward explicit schemes for 
discretization. However, while in the two scalar mod- 
els (SM, RRV) about the same maximum time step is 
possible, the tensorial models allow larger time steps; a 
simulation with the LCT model gains a factor of about 
ten in time steps over the scalar models. While by use 
of adaptive mesh techniques 24J (and implicit schemes), 
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the overall running time can certainly be reduced by more 
than this factor for large systems, the advantage of the 
tcnsorial models may persist even in such a setting as it 
is consistently present in a range of grid spacings. 

Next, it is interesting to compare how the different 
models behave regarding their relaxation to a stationary 
profile when initialized with a straight interface having a 
width that is either too small or too large. These simu- 
lations are done without elasticity, i.e., for F — 0. First, 
we verify that all the models remain in their equilibrium 
state when initialized with a tanh profile of the correct 
width. 
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FIG. 6: Relaxation towards the correct interface profile for 
the SM model. The initial interface width is 0.25 W, the time 
/ is given in units of W 4 /M. 

For brevity of language, we define here a profile 
tanhz/W to have width W, even when the width on which 
it rises from -0.9 to 0.9 rather is 2.94 W. In all plots, the 
time / is given in units of W 4 /M, and only a small section 
about the interface is shown. 

All of the models do reasonably well in relaxing from 
a planar profile of width 0.25 W to their equilibrium 
state, see Figs. [6] to [9j In our implementation and with 
the given sets of parameters, simulations with the RRV 
model broke down if the initial interface width was cho- 
sen to be smaller than 0.23 W. The RRV and GCT model 
relax quickly to a final profile of width W, while the SM 
model needs a little more time (but the permissible time 
step is larger for the GCT model, so the numerical runnig 
time is shortest for it). In addition, some care has to be 
taken in the discretization to make the LCT model deal 
efficiently with too thin interfaces. 

To see this, we write Q 4 = 1 - n : n + bin : n with 



w 2 (v<py 

4/(0 



(107) 




FIG. 7: Relaxation towards the correct interface profile for 
the RRV model. The initial interface width is 0.25 W, the 
time t is given in units of W 4 /M. 



5 (1 — tanhz/£) we find 



(108) 



Taking an interface of width £ with the profile cf>(z) — 



which becomes very large for £ <K W. In fact, for 
£, = 0.25VK, we have b A Q = 50625, a number that, when 
plugged into the equations of motion, would impose a 
prohibitively small time step for stability (or accuracy, 
in an implicit scheme). Hence, we introduce a cutoff 
for on the order of 50. In production runs, where 
one normally starts with a front profile having at least 
approximately the correct width, a cutoff of 10 may be 
sufficient. 

When the profile is initialized with too large a width 
[Figs. 10 through |13] , more interesting differences can 
be seen. Not unexpectedly, the GCT model [Fig. 13 
is the one making the least fuss about an interface five 
times too wide: that the model is nonconservative on the 
scale where the phase field varies strongly is an advantage 
here. The interface approaches its correct width in a 
time of about t « 1.25 W 2 /N, which corresponds to f ~ 
1 W 4 /M for our parameter choice. For the other models, 
this takes much longer, as this kind of adaptation requires 
diffusion orthogonally to the front, which is slow because 
it is suppressed in the asymptotic limit. 

The RRV model goes through a series of transforma- 
tions of the profile involving as an intermediate state a 
spatially varying slope in the vicinity of the contour line 
0=1 defining the interface position. Even after a time 
off* 30W 4 /M, while near the interface position the pro- 
file is well-behaved and has the right width, there are 
still indentations in it far from the interface, and these 
disappear only slowly. 

While the LCT model keeps a nicer profile all the time, 
it relaxes only slowly as well. Moreover, if the boundary 
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FIG. 8: Relaxation towards the correct interface profile for 
the LCT model. The initial interface width is 0.25 W, the 
time t is given in units of W 4 /M. 



FIG. 10: Relaxation towards the correct interface profile for 
the SM model. The initial interface width is 5 W, the time t 
is given in units of W 4 /M. 
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FIG. 9: Relaxation towards the correct interface profile for 
the GCT model. The initial interface width is 0.25 W, N = 
1.25 M/W 2 , the time t is given in units of W 4 /M. 



FIG. 11: Relaxation towards the correct interface profile for 
the RRV model. The initial interface width is 5 W, the time t 
is given in units of W 4 /M. 



values of cfi are not fixed to be equal to zero or one, it will 
relax to constant values in the bulk different from these 
ideal values (in the absence of elasticity) . Indeed, inspec- 
tion of Eq. (64 1 shows immediately that any constant 
value of (f> solves the bulk equations of motion. (This is 
no longer true in the presence of elasticity.) For phases 
extending to the system boundary, the value of the con- 
stant is only fixed by the boundary conditions. Therefore, 
the model should always be run with Dirichlct boundary 
conditions for the phase field. (Due to the conservation 
law, inclusions of one phase in another will keep their 
value, even in the presence of elasticity, if correctly 
initialized to zero or one, as long as their inner volume 
is much larger than that of their interface.) Performing 



such a simulation, we found relaxation to be as slow as 
for the SM and RRV models but the interface profile to 
look more reasonable. 

To summarize, when interface thickness is believed to 
be an issue in simulations, i.e., when there are reasons 
to think that it might vary considerably (which may be 
the case when surface tension anisotropy is included in 
the model), the nonexact realization of the conservation 
law by the GCT model may turn out a virtue rather 
than a drawback, since changes in the direction normal 
to the interface by diffusion only, as realized in the other 
models, tend to be too slow. 

Finally, we compare the different dynamics for a "real- 
life" situation of an elliptical inclusion that morphs into 
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tM/W 4 



FIG. 12: Relaxation towards the correct interface profile for 
the LCT model. The initial interface width is 5 W, the time t 
is given in units of W 4 /M. 



FIG. 14: Comparison of the time development of the size 
of an elliptical inclusion. The square system had the length 
L/W = 20. The initial ellipse had a semimajor of flo = L/4 
and a semiminor of L/8. All models except the SM-model 
converge to circles with the same radius r. 




FIG. 13: Relaxation towards the correct interface profile for 
the GCT model. The initial interface width is 5 W, the time 
/ is given in units of W 4 /M. 



a circle without elastic effects, F — 0. The system is 
initialized with a sharp interface ellipse with semimajor 
ao and semiminor «o/2 and is then allowed to relax for 
a few time steps running the GCT-model (with a La- 
grange multiplier A — 0), in order to obtain an initial 
condition with the correct interface width everywhere. 
We then measure the time evolution of the semimajor 
and semiminor of the ellipse, continuing the run with the 
model to be studied. 

As Fig. [14] shows, all models but the SM-model con- 
verge to a circle with the correct radius V2«o/2. The SM- 
model shows different behavior, namely a too small ra- 
dius that seems to decrease further. Since the phase field 
is a conserved quantity (we also checked that numerically 



for our code), this can only mean that the final shape of 
the inclusion is not a true circle but a slightly deformed 
one, displaying a certain level of anisotropy. We then in- 
crease the size of the system and the included ellipse while 
keeping the interface width constant, resulting in a better 
scale separation oq/W. While for the LCT-, GCT- and 
RRV- model the curves collapse onto a single line, this 
is not the case for the SM-model. Fig. [15] demonstrates 
this behavior for the LCT-model. The comparison for 
the SM-model is shown in Fig. [16] 
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FIG. 15: Elliptical inclusion: comparison of the time devel- 
opment of the length of the semimajor a for the LCT-model. 
The initial length is denoted by ciq and the different curves 
correspond to different scale separations a /W. All curves 
collapse onto a single line. 
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FIG. 16: Elliptical inclusion: comparison of the time devel- 
opment of the length of the semimajor a. The initial length is 
denoted by a and the different curves correspond to different 
scale separations ag/W. The performance of the SM-model 
becomes asymptotically better for larger systems. Results for 
the LCT-model have been included as a reference. 



VIII. CONCLUSIONS 

The intuitive approach to constructing a phase-field 
model for surface diffusion consists in using the chemical 
potential known from the nonconservative model to de- 
fine a current, involving its gradient and a mobility that 
vanishes in the bulk phases, and in taking the divergence 
of this current as the time derivative of the phase field. 
As has been shown in this article, this approach - the 
SM model - fails to produce the correct asymptotics in 
a subtle way. It does reproduce the equlibrium limit cor- 
rectly and it appears to work numerically, although less 
efficiently than the alternatives discussed. 

We offer a simple argument why the SM model should 
not be expected to work properly: The chemical potential 
functional of the model is constructed so that the chem- 
ical potential vanishes in the bulk phases. As the dif- 
fusion operator is essentially a scalar, diffusion acts also 
orthogonally to the interface; in its vicinity, the effect is 
even strong, because the slope of the phase field is largest 
in the direction perpendicular to the corresponding level 
set. This diffusive effect constitutes a driving force for re- 
laxation of the chemical potential towards zero also close 
to the interface (asymptotically, the chemical potential 
is zero at next-to leading order) . Surface diffusion of the 
chemical potential is then not the only effect contributing 
to the interface dynamics. 

The RRV model avoids this problem by leaving the 
chemical potential in the bulk undetermined. Absence of 
diffusion in the bulk is not guaranteed by the chemical 
potential but by the vanishing mobility. Since the bulk 
chemical potential is free to vary, a true interface chemi- 
cal potential can build up, the surface diffusion of which 
governs the interface dynamics. 



Our contribution in this article is to explore the idea 
that the failure of the SM model might be remedied in- 
stead by making the mobility a tensor. After all, sur- 
face diffusion may be interpreted as highly anisotropic 
three-dimensional diffusion with a diffusion tensor that 
has zero eigenvalue in one direction. Whereas the pre- 
existing RRV model has the correct asymptotics, it does 
not exploit that idea. A straightforward attempt of its 
realization however fails in a rather drastic way, because 
restricting diffusion to the surfaces of constant phase field 
does not impose any functional dependence of cf> in the 
normal direction given by this foliation. 

Modifying the tensorial mobility, one obtains the LCT 
and GCT models, both exhibiting the correct asymptotic 
behavior. 

An analytic linear stability analysis of a planar front 
demonstrates that these two models reproduce the cor- 
rect dispersion relation for local perturbations of the pro- 
file corresponding to a modified interface position. The 
phase-field dispersion relation is even free of corrections 
due to the front width. A similar analysis turns out in- 
feasible for the SM model, while it could in principle be 
completed numerically for the RRV model (after analytic 
reduction to an ordinary differential equation). 

Numerical study of the four models suggests that 
whereas the SM model has a range of quantitative va- 
lidity (despite its not being asymptotic), and the RRV 
model is definitely viable, the modified tensorial models 
discussed here are probably more useful for large-scale 
simulations, as they permit larger time steps. 

A possible way to understand this is as follows. In the 
scalar-mobility models (SM and RRV), the inner equa- 
tion determining the interface velocity appears only three 
orders in W after the leading order. A simulation must 
of course represent the full model equations. Suppose we 
wish the interface velocity to be determined with an error 
not exceeding order W. Then the leading-order equation 
must be simulated with an accuracy 0(W 4 ) at least. If 
one takes a simple second-order accurate discretization 
for the gradient operators in the equations, the numeri- 
cal error due to discretization alone will be 0(Ax 2 ) for a 
grid spacing Ax. To keep this smaller than W 4 , Ax would 
have to scale with W 2 , which can be expected to lead to 
large computation times. Therefore, in the scalar mobil- 
ity models, high-accuracy discretizations are mandatory, 
even when the asymptotic error is controllable, as is the 
case for RRV. 

On the other hand, in the LCT and GCT models, dis- 
VAlandlVB 



cussed in Sees. 



as soon as the phase-field 



profile is represented with an error of order W or better, 
the effective leading order is only one order lower than 
the one determining the interface velocity, similar to the 
nonconservative case. Hence, reasonable accuracy should 
be attainable with grid spacings that scale as W, not as 
W 2 . None of the benefits of high-accuracy discretizations 
would be lost to the need of representing terms very accu- 
rately that are very small in the leading-order equation. 
Generally speaking, the most efficient model in terms 
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of computational cost is the LCT model. However, it is 
slowest in terms of "real" time (though not in terms of 
the number of time steps), when it comes to the relax- 
ation from a wrong width of a planar interface to the 
correct value. This may be seen as a signature that the 
model is most efficient in suppressing diffusion perpen- 
dicular to the interface. In optimally initialized sim- 
ulations, interface width variations arise gradually, via 
curvature changes and/or orientation changes (in models 
with anisotropic surface tension) and one should expect 
their relaxation via diffusion along the interface to be 
sufficiently efficient. Nevertheless, this property of the 
LCT model reduces its robustness as a numerical tool. 

This is why we suggest the GCT model as an alterna- 
tive that seems to be more accurate in most applications 
than the LCT model and digests variations of the inter- 
face width more easily than all the other models. Both 
accuracy and robustness of the GCT model, connected 
with its still favorable efficiency and simplicity of imple- 
mentation, should make it the approach of choice in most 



where the derivatives are to be taken for r — > +0, should 
they be discontinuous at r — 0. Analogous expressions 
with r — > -0 are obtained for the asymptotics as p — > -oo. 

Equating powers of W, we then successively get the 
asymptotic relationships 



lim Y (p) = «Ao(±0) , 

p— >±oo 

¥i(p) ~ P«Ao(±0) + <Ai(±0) (p -» ±oo) , 



¥ 2 (p) 



ipV '(±0)+p^(±0) + ^ 2 (±0) 



(P 



3). 



(A4) 



(A5) 



(A6) 



Moreover, asymptotic relations such as (A5l can be de 



composed into statements about function limits 
lim 5 p Ti(p) = ^(±0), 

p— >±oo 

lim K(p)-p^(±0)l = 0i(±O). 

p— >±oo L J 



(A7) 
(A8) 
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APPENDIX A: MATCHING CONDITIONS 

Let ij/(x, z, i) — if/(r, s, t) be some arbitrary (sufficiently 
often differentiable) function of space and time obtained 
in solving the outer equations. We write the correspond- 
ing function of the inner solution as ^(p, s, ?), and sup- 
press from now on, in this section, the dependence of 
functions on s and t. Moreover, we write the coefficient 
functions in expansions with respect to W with simple 
subscripts indicating their order rather than superscripts 
in parentheses as in the main text. There, the notation 
is dictated by the fact that a subscript would interfere 
with certain other subscripts; here, a superscript would 
interfere with the primes denoting derivatives. 

We must have the asymptotic relationship 

¥(p) ~ tAO) = t/r(Wp) (p -> oo, W -> 0, Wp -> 0) . (Al) 

Expanding both functions in powers of W, we get 

¥(p) = <F (p) + (p) + W 2 y 2 (p) + ..., (A2) 
0(Wp) = Mr) + Wiffi(r) + W 2 iff 2 (r) + ... 

= M0) + wiM(0) + iM°)] 

+ wV^o (0) + P^(0) + <A 2 (0)] + . . . , 



APPENDIX B: USEFUL PROPERTIES OF THE 
PHASE FIELD FUNCTIONS 

In order to simplify it for the reader to find the ac- 
tual relationships for the various functions involving the 
phase-field that are used in the text, they are collected 
here for reference (and concreteness). Often only certain 
properties but not the precise form of these functions are 
important. 

The chosen double-well potential is 



/(0) = 2 (l-0) 2 . 
Its derivative is given by 

/'(<*) = 20(1 - 0)(1 - 2<f) , 
its second derivative reads 

f"{<j>)= 2 [1-60(1-0)] 



(Bl) 



(B2) 



(B3) 



Let us define another function h(<f>), which turns out use- 
ful, by 



/i(0) = 2 (3-20), 
having the derivative 

= 60(1-0). 



(B4) 



(B5) 



To solve the ordinary differential equation satisfied by 
the zeroth-order inner solution 



<3 PP <D (0) - 2/'(O (0) ) = , 



(B6) 



(A3) 



with boundary conditions limp^-a, <t> (0) (p) = 1 and 
linip^oo 4> (0) (p) = 0, we multiply by (9 P <E> (0) , integrate and 
take the square root (with the correct sign) to obtain 

<9p(D (0) = -20> (0, (1 - <D (0) ) = -i//(<£ (0) ) , (B7) 
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which can be solved by separation of variables. The so- 
lution is, up to a translation in p, given by 



O (0) = -(1 -tanhp) . 



(B8) 



Requiring the position of the interface to be at p — 
fixes the choice out of the one-parameter set of solutions, 
present due to the translational invariance of the differ- 
ential equation. 

With the help of the second equality of ( |B7[ ) , it is easy 
to calculate certain integrals appearing in the asymptotic 
analysis. Those integrals typically contain the factor 

(<9 p <t> (0) ) ; to do the integral, it is then beneficial to re- 
place one of the factors (and only one) with -/z'(cD (0) )/3. 
Integrals obtained this way have the structure 

CO 

1 r°° 

= -- dp/(/l(O (0) ))/z'(O (0) )^ p <^> (0, 

■J %J —CO 

3 Ji 

i r 1 

= - dhofih,). (B9) 
J Jo 

This way, one arrives, for example, at 
' i 1 

a,»") 2 * = i. 



J 



(BIO) 



In the RRV model [23], we have to evaluate two more 
integrals, namely 



g(d> {0) ) d p <fr (0) dp = - ( 10x 2 (l - xfdx = -1 

oo Jo 

Xoo pea 
B(®<®)dp= 3(3 P <D (0) ) =1. 
OO v7 — CO 



/3. 



(Bll) 



we then have, at leading order 



W-2^o = -4Mi// 



(C4) 



meaning that cither w_2 = -4M or i^o = 0. Thus we have 
to distinguish two cases. 

If w_2 - — 4M, the lowest-order inner equation becomes 



TO • 



cosh 2 Z 



^o(Z) = , 



(C5) 



where we have used ( |B8| . This is a one-dimensional 
Schrodinger equation that can be reduced to Legendre's 
differential equation [32] via the substitution u = tanhZ 
and hence is exactly solvable. In this particular case, the 
solution that remains finite at infinity is the second-order 
(in u) Legendre polynomial 



¥ (Z)= ^(3tanh 2 Z-l) . 



(C6) 



This means that the perturbation is not localized - it does 
not approach zero for Z — > +oo. As usual, once the zeroth- 
order solution has been found, perturbation theory can 
be carried through without major difficulties, inserting 
the expansions ( C3 1 into (CI) and ( C2 1 . Calculating the 



eigenvalue to zeroth order in W, we find a> a from ( 83 1 



In fact, continuing the calculation to higher orders, we 
note that no additional contributions to a> a arise. The 



suspicion that ( C6 ) is an exact solution - and u> a the 



corresponding exact eigenvalue to Eq. (CI I is confirmed 
by backsubstitution into the equation. 

The second solution is obtained by assuming «_2 = 
and, hence, i^o(z) — in the outer region. This gives the 
leading-order inner equation 



APPENDIX C: DETAILS OF THE ANALYTIC 
LINEAR STABILITY CALCULATION 

In the nonconservative case, the eigenvalue problem to 
be solved reads 

M 

This may be considered the inner problem, given in the 
coordinates (x,Z), with Z = z/W. Setting if/(z) = ^(Z), the 
corresponding outer problem reads: 



^ = §2 i dzz ~ W2kl ~ 2 /"(°o)) ¥ 



(CI) 



4 \ 



(C2) 



*F '(Z)-2/"(a>oW>(Z) = 0, 



(C7) 



the relevant solution of which we know already, because 



the linear operator acting here is just £. from (39 1 with 
the role of p taken by Z. Hence 



¥ (Z) = %{Z) = - 



1 



2 cosh 2 Z 



(C8) 



which is a solution localized about the interface and ap- 
proaching zero for Z — > +oo. Inserting the expansions 
( C3 1 with the evaluated results for oj-2, ^o(z) and 'I'o(Z) 



where we have used lim z ^ ±M /"(O ) = 2 [see |B3] |]. Using into |C1| and (|C2j, we find that the eigenvalue is to all 
the expansions 

CO-2 W_i 



orders given by cot, from (83 1 and that, in fact, ^oiZ) and 



T(Z) = W (Z) + W¥i(Z) + W 2x ¥ 2 (Z) + . . 
Hz) = Mz) + Wi(z) + wV 2 (z) + . . . 



(C3) 



a>b constitute an exact solution to the eigenvalue prob- 
lem. There may be more admissible solutions, but for all 
of them a) should be negative and diverge more strongly 
than 1 /W 2 for W — > 0, so they will not be relevant in the 
limit of small W. 
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